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f^ Abstract. We present the new general-relativistic magnetohydrodynamics (GRMHD) 

^^ capabilities of the Einstein Toolkit, an open-source community-driven numerical rel- 

. . ativity and computational relativistic astrophysics code. The GRMHD extension of 

^ J^ the Toolkit builds upon previous releases and implements the evolution of relativis- 

^^ tic magnetised fluids in the ideal MHD limit in fully dynamical spacetimes using the 

Jh same shock-capturing techniques previously applied to hydrodynamical evolution. In 

order to maintain the divergence-free character of the magnetic field, the code imple- 
ments both hyperbolic divergence cleaning and constrained transport schemes. We 
present test results for a number of MHD tests in Minkowski and curved spacetimes. 
Minkowski tests include aligned and oblique planar shocks, cylindrical explosions, mag- 
netic rotors, Alfven waves and advected loops, as well as a set of tests designed to study 
the response of the divergence cleaning scheme to numerically generated monopoles. 
We study the code's performance in curved spacetimes with spherical accretion onto 
a black hole on a fixed background spacetime and in fully dynamical spacetimes by 
evolutions of a magnetised polytropic neutron star and of the collapse of a magnetised 
stellar core. Our results agree well with exact solutions where these are available and 
we demonstrate convergence. All code and input files used to generate the results are 
available on http://einsteintoolkit.org. This makes our work fully reproducible 
and provides new users with an introduction to applications of the code. 

PACS numbers: 04.25.D-, 04.30.-w, 04.70.-s, 07.05. Tp, 95.75.Pq 
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1. Introduction 

The recent years have seen rapid developments in the field of numerical relativity. 
Beginning with the first fully general-relativistic (GR) simulations of binary neutron 
star (NS) mergers by Shibata and Uryu in 1999 pQ, fully dynamical general-relativistic 
hydrodynamics (GRHD) has been explored by a growing number of research groups. 
A major step forward occurred 2005, when three independent groups developed 
two different techniques, the generalised harmonic gauge formalism [2] and the so- 
called "moving puncture" method [21 S], to evolve vacuum spacetimes containing 
black holes (BH) without encountering numerical instabilities that mark the end of 
a numerical simulation. Since then, these evolution schemes have been incorporated 
in hydrodynamic simulations and now one or the other is used in essentially every full 
GR calculation of merging binaries, as well as other astrophysical phenomena involving 
dynamical spacetimes, such as stellar collapse and BH formation (e., [5]). 

A more recent advance concerns the incorporation of magnetic fields, particularly 
in the ideal magnetohydrodynamics (MHD) limit, into dynamical simulations. Building 
upon work that originally focused on astrophysical systems in either Newtonian gravity 
or fixed GR backgrounds, as is appropriate for studies of accretion disks (see [S]- 
[H] for reviews of the topic), MHD evolution techniques have been incorporated into 
dynamical GR codes and used to study the evolution of NS-NS (e.g., P-IT^) and BH- 
NS mergers (e.g., [131 13)) self-gravitating tori around black holes [ISIIIS], neutron star 
collapse (e.g., [I71[TB]), stellar core collapse [IH], and the evolution of magnetised plasma 
around merging binary BHs (e.g., [201 EI])- Given their rapidly maturing capabilities 
and widespread use, many numerical relativity codes have been made public, either in 
total or in part. The largest community-based effort to do so is the Einstein Toolkit 
consortium [221 [23], which has as one of its goals to provide a free, publicly available, 
community-driven general-relativistic MHD (GRMHD) code capable of performing 
simulations that include realistic treatments of matter, electromagnetic fields, and 
gravity. The code is built upon several open-source components that are widely used 
throughout the numerical relativity community, including the Cactus computational 
infrastructure [21H2S], the Carpet adaptive mesh refinement (AMR) code [23"E2], the 
McLachlan GR evolution code [SUl EI], and the GRHD code GRHydro [22j, which has 
been developed starting from a public variant of the Whisky code [32H3S]- Many of the 
details of the Einstein Toolkit may be found in [22j, which describes the routines used to 
provide the supporting computational infrastructure for grid setup and parallelization, 
constructing initial data, evolving dynamical GRHD configurations, and analysing 
simulation outputs. An extension of the Einstein Toolkit to general multi-patch grids 
with Cartesian and curvilinear geometry is discussed in [30]- 

In this paper, we describe the newest component of the Einstein Toolkit: a MHD 
evolution scheme incorporated into the GRHydro module. Using what is known as the 
Valencia formulation [S7HS5] of the GRMHD equations, the new code evolves magnetic 
fields in fully dynamical GR spacetimes under the assumptions of ideal MHD, i.e., the 
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resistivity is taken to be zero and electric fields vanish in the comoving frame of the fluid. 
While the evolution equations themselves are easily cast into the same flux- conservative 
form as the other hydrodynamics equations, the primary challenge for evolving magnetic 
fields is numerically maintaining the divergence-free constraint of the magnetic field. For 
this, we have implemented a variant of the constrained transport approach |1D] and a 
hyperbolic "divergence cleaning" technique, similar to that discussed in [¥T| H2]. 

The GRHydro module of the Einstein Toolkit is a new, independently written 
code, and — in particular — the development is completely independent of the WhiskyMHD 
code [43j, since the GRHydro development started from the GRHD (non-MHD) version 
of the Whisky code. Given its close interaction with the hydrodynamic evolution, 
the MHD routines are contained within the GRHydro package (or "thorn", in the 
language of Cactus-based routines). It is developed in a public repository, along 
with other components of the Einstein Toolkit, and it is publicly available under the 
same open-source licensing terms as other components in the Toolkit. Given the 
public nature of the project, the code release includes the subroutines themselves, 
complete documentation for their use, and parameter files needed to reproduce the 
tests shown here. Code and parameter files are available on the Einstein Toolkit web 



page, http : //einsteintoolkit . org 



In order to validate our GRMHD implementation, we include a number of flat- 
space tests whose solutions are either known exactly or approximately, and which 
have frequently been used as a testbed for other GRMHD codes. These include 
planar and cylindrical shocktubes, a rotating bar threaded by a magnetic field (the 
"magnetic rotor"), propagating Alfven waves, and the advection of a flux loop. We 
test the code's ability to handle static curved spacetime by studying Bondi accretion 
onto a Schwarzschild BH and perform GRMHD simulations with dynamical spacetime 
evolution for a magnetised TOV star and for the collapse of a rotating magnetised stellar 
core to a rotating neutron star. 

This paper is structured as follows: in section |2} we describe the Einstein Toolkit 
and briefly discuss other relativistic MHD codes. In section |3j we describe the Valencia 
formulation of the GRMHD equations, and in section |4] the numerical techniques used 
in GRHydro. In section [5| we describe the tests we have carried out with the new code. 
Finally, in section [6} we summarise and discuss future directions of the Einstein Toolkit. 

2. Numerical Relativistic Hydrodynamics and GRMHD codes 

The GRHydro code is, to the best of our knowledge, the first publicly released 3D 
GRMHD code capable of evolving configurations in fully dynamical spacetimes. Still, 
its development makes use of code and techniques from a number of other public codes, 
which we summarise briefly here. The most obvious of these is the Einstein Toolkit, 
within which its development has taken place. As we discuss below, the Toolkit itself 
is composed of tens of different modules developed by a diverse set of authors whose 
number is approaching 100. The MHD techniques, while independently implemented in 
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nearly all cases, rely heavily on both the numerical techniques found within the GRHydro 
package [22], and thus the original Whisky code [32HSS], as well as the publicly available 
HARM code [il] . 

2.1. The Einstein Toolkit 

The Einstein Toolkit includes components to aid in performing relativistic astrophysical 
simulations that range from physics modelling (initial conditions, evolution, analysis) 
to infrastructure modules (grid setup, parallelization, I/O) and include related tools 
(workflow management, file converters, etc.). The overall goal is to provide a set of 
well-documented, well-tested, state-of-the-art components for these tasks, while allowing 
users to replace these components with their own and/or add additional ones. 

Many components of the Einstein Toolkit use the Cactus Computational 
Toolkit [21H2n], a software framework for high-performance computing. Cactus 
simplifies designing codes in a modular ("component-based") manner, and many existing 
Cactus modules provide infrastructure facilities or basic numerical algorithms such as 
coordinates, boundary conditions, interpolators, reduction operators, or efficient I/O in 
different data formats. 

The utilities contained in Einstein Toolkit, e.g., help manage components [151 HB] . 
build code and submit simulations on supercomputers [HJ HH], or provide remote 
debuggers [IH] and post-processing and visualisation interfaces for Visit [50^ . 

Adaptive mesh refinement (AMR) and multi-block methods are implemented by 
Carpet [27H2U] which also provides MPI parallelization. Carpet supports Berger-Oliger 
style ("block-structured") adaptive mesh refinement [SI] with sub-cycling in time as 
well as certain additional features commonly used in numerical relativity (see [23 fo'^ 
details). Carpet provides also the respective prolongation (interpolation) and restriction 
(projection) operators to move data between coarse and fine grids. Carpet has been 
demonstrated to scale efficiently up to several thousand cores [22j- 

Carpet supports both vertex-centred and cell-centred AMR. In vertex-centred 
AMR, each coarse grid point (cell Center) coincided with a fine grid point (cell Center), 
which simplifies certain operations (e.g. restriction, and also visualisation). In this 
paper, we present results obtained with vertex-centred AMR only. In cell-centred 
AMR, coarse grid cell faces are aligned with fine grid faces, allowing in particular exact 
conservation across AMR interfaces. For a more detailed discussion of cell-centred AMR 
in hydrodynamics simulations with the Einstein Toolkit we refer the reader to [5B] . 

The evolution of the spacetime metric in the Einstein Toolkit is handled by the 
McLachlan package [201 EI]- This code is autogenerated by Mathematica using the 
Kranc package [S2HS1], implementing the Einstein equations via a 3 + 1-dimensional 
split using the BSSN formalism [55145^ . The BSSN equations are finite-differenced at a 
user-specified order of accuracy, and coupling to hydrodynamic variables is included via 
the stress-energy tensor. The time integration and coupling with curvature are carried 
out with the Method of Lines (MoL) [20], implemented in the MoL package. 
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Hydrodynamic evolution techniques are provided in the Einstein Toolkit by the 
GRHydro package, a code derived from the public Whisky GRHD code [32HS5]- The 
code is designed to be modular, interacting with the vacuum metric evolution only 
by contributions to the stress-energy tensor and by the local values of the metric 
components and extrinsic curvature. It uses a high-resolution shock capturing finite- 
volume scheme to evolve hydrodynamic quantities, with several different methods for 
reconstruction of states on cell interfaces and Riemann solvers. It also assumes an 
atmosphere to handle low-density and vacuum regions. In particular, a density fioor 
prevents numerical errors from developing near the edge of matter configurations. 
Boundaries and symmetries are handled by registering hydrodynamic variables with 
the appropriate Cactus routines. Passive "tracers" or local scalars that are advected 
with the fluid flow may be used, though their behaviour is generally unaffected by the 
presence or absence of magnetic fields. 

In adding MHD to the pre-existing hydrodynamics code, our design philosophy has 
been to add the new functionality in such a way that the original code would run without 
alteration should MHD not be required. For nearly all of the more complicated routines 
in the code, this involved creating parallel MHD and non-MHD routines, while for some 
of the more basic routines we simply branch between different code sections depending 
on whether MHD is required or not. While this requires care when the code is updated, 
since some changes may need to be implemented twice, it does help to insure that users 
performing non-MHD simulations will be protected against possible errors introduced 
in the more frequently changing MHD routines. 

2.2. Other relativistic MHD codes 

In extending the Einstein Toolkit to include MHD functionality, we incorporated 
techniques that have previously appeared in the literature. Of particular importance is 
the HARM code [lH EB E2], a free, publicly available GRMHD code that can operate on 
fixed background metrics, particularly those describing spherically symmetric or rotating 
black holes. Our routines for converting conservative variables back into primitive ones 



are adapted directly from HARM ([63]; see section 4.4 below), as is the approximate 
technique used to calculate wave speeds for the Riemann solver ([H]; see section 4.3 
below) . 

Since many groups have introduced relativistic MHD codes, there are a number 
of standard tests that may be used to gauge the performance of a particular code. 
Many of these are compiled in the descriptions of the HARM code mentioned above, as 
well as papers describing the Athena MHD code [Ml ES] , the Echo GRMHD code [HS] , 
the Tokyo/Kyoto group's GRMHD code [67], the UIUC GRMHD code §SiW\, the 
WhiskyMHD code |13], and the LSU GRMHD code [70]. We have implemented several of 
these tests, covering all aspects of our code, as we describe in section |5] below. 
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3. The Valencia formulation of ideal MHD 

While GRHydro does not technically assume a particular evolution scheme to evolve 
the field equations for the GR metric, the most widely used choice for the ET code 
is the BSSN formalism [501 EZ], the particular implementation of which is provided by 
the McLachlan code [301 EI]- We do assume that the metric is known in the ADM 
(Arnowitt-Deser-Misner) form [71J, in which we have 

ds^ = g^^dxf'dx'' = (-a^ + (3i/3')dt'^ + 2(3idt dx' + -fijdx'dx^ , (1) 

with g^^, a, (5\ and 7ij being the spacetime 4-metric, lapse function, shift vector, and 
spatial 3-metric, respectively. Note that we are assuming spacelike signature, so the 
Minkowski metric in fiat space reads rj^^ = diag(— 1, 1, 1, 1). Roman indices are used for 
3-quantities and 4-quantities are indexed with Greek characters. We work in units of 
c = G = Mq = 1 unless explicitly stated otherwise. 

GRHydro employs the ideal MHD approximation - fiuids have infinite conductivity 
and there is no charge separation. Thus, electric fields E^, = u^E'^^'^ in the rest frame 
of the fiuid vanishes and ideal MHD corresponds to imposing the following conditions: 

u.E^-' = . (2) 

Note that throughout this paper we rescale the magnitude of the relativistic Faraday 
tensor E^'^ and its dual y^^" = ^e^'^'^'^E^x, as well as the magnetic and electric fields, by 
a factor of 1/^/4tt to eliminate the need to include the permittivity and permeability of 
free space in cgs-Gaussian units. 

Magnetic fields enter into the equations of hydrodynamics by their contributions 
to the stress-energy tensor and through the solutions of Maxwell's equations. The 
hydrodynamic and electromagnetic contributions to the stress-energy tensor are given, 
respectively, by 

T^'' = phu^'u" + Pg^" = {p + pe + P) u^v!" + Pg^"" , (3) 

and 

T^^ = E^^E'^x - \9^''F^-Ex. = b'^u^u'^ - h^lf + ^-g^^ , (4) 

where p, e, P, u^ ^ and h = 1 + e + P j p are the fiuid rest mass density, specific 
internal energy, gas pressure, 4-velocity, and specific enthalpy, respectively, and h^ is 
the magnetic 4- vector (the projected component of the Maxwell tensor parallel to the 
4-velocity of the fiuid): 

6^ = n/F^' . (5) 

Note that l? = U^h^ = 2Pm, where Pm is the magnetic pressure. Combined, the stress- 
energy tensor takes the form: 

T'"' = (^p + pe + P + h^) u^v!' + ( ^ + Y ) ^''^ ~ ^''^^ (^) 
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where we define the magnetically modified pressure and enthalpy, P* = P + Pm = 
P + 572 and /I* = 1 + e + (P + h"^) /p, respectively. 

The spatial magnetic field (living on the spacelike 3-hypersurfaces), is defined as 
the Eulerian component of the Maxwell tensor 

B' = n^p'^ = -aY'' , (7) 

where n^ = [—a, 0, 0, 0] is the normal vector to the hypersurface. 

The equations of ideal GRMHD evolved by GRHydro are derived from the local GR 
conservation laws of mass and energy-momentum, 

V^.r = , V^T^" = , (8) 

where V^.^ denotes the covariant derivative with respect to the 4- metric and J^ = pu^ 
is the mass current, and from Maxwell's equations, 

v;f'' = . (9) 

The GRHydro scheme is written in a first-order hyperbolic fiux-conservative 
evolution system for the conserved variables D, S^, r, and i3*, defined in terms of the 
primitive variables p,e,v\ and 5* such that 

D = ^pW , (10) 

S^ =^(^ph*W\j - ab%j) , (11) 

r = ^ (^ph*W^ - P* - (aby) - D , (12) 

B' = VlB' , (13) 

where 7 is the determinant of jij. We choose a definition of the 3- velocity f* that 
corresponds to the velocity seen by an Eulerian observer at rest in the current spatial 
3-hypersurface [72], 

-^-^ 

and W = {1 — v^Vi) ^1"^ is the Lorentz factor. Note that f*, B'\ S\ and /3* are 3- vectors, 
and their indices are raised and lowered with the 3-metric, e.g., Vi = 'jijV^ . 

The evolution scheme used in GRHydro is often referred to as the Valencia 
formulation [HI I57H5U] . Our notation here most closely follows that found in [75] . 
The evolution system for the conserved variables, representing ([s]) and the spatial 
components of ([O]), is 



with 



U = [D,S„T,B>'], 

Dv' 

Sjv' + y/^P*6;-bjByW 

TV' + J^P*v' - ab^B'/W 

B^v' — B'v'^ 



a X 



(16) 
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OCy/l -X 






0, 



8 In a T-i/ij/pO 





(17) 



Here, w* = f * — /3Yq; and V^^^ are the 4-Christoffel symbols. 

The time component of ([9]) yields the condition that the magnetic field is divergence- 
free, the "no-monopoles" constraint: 

^ dA^B^)=Q, (18) 



V -B 



^ 



which also implies 

diB' = . 



(19) 



In practice, we implement two different methods to actively enforce this constraint. In 
the "divergence cleaning" technique we include the ability to modify the magnetic field 
evolution by introducing a new field variable that dissipates away numerical divergences, 
which we discuss in detail in section 4.5. 1| below. An alternative method, commonly 
called "constrained transport", instead carefully constructs a numerical method so that 



the constraint (19) is satisfied to round-off error at the discrete level. We discuss this 



method in section 14.5.21 



4. Numerical Methods 

GRHydro 's GRMHD code uses the same infrastructure and backend routines as its 
pure general relativistic hydrodynamics variant |22]- In the following, we focus on the 
discussion of the numerical methods used in the extension to GRMHD. 

GRHydro's GRMHD code implements reconstruction of fluid and magnetic fleld 
variables to cell interfaces for all of the methods present in the original GRHydro code: 
TVD (total variation diminishing) (e.g., [Zl]), PPM (piecewise parabolic method) [75], 
and ENO (essentially non-oscillatory) [751 177]. In addition, we added enhanced PPM 
(ePPM, as described in [23 EH]), WEN05 (5*^^ order weighted-ENO) |72], and MP5 
(5*^ order monotonicity preserving) [SD] both in GRMHD and pure GR hydrodynamics 
simulations. We discuss the different reconstruction methods in section 14.21 The 
code computes the solution of the local Riemann problems at cell interfaces using the 
HLLE (Harten-Lax-van Leer-Einfeldt) [HH IH2] approximate Riemann solver discussed 
in section |4.3[ We implement the conversion of conserved to primitive variables for 
arbitrary equations of state (EOS), including polytropic, F-law, hybrid polytropic/P- 
law, and microphysical flnite-temperature EOS. We summarise the new methods for 



GRMHD in section 4.4 An important aspect of any (GR)MHD scheme is the numerical 



method used to preserve the divergence free constraint. GRHydro implements both the 
hyperbolic divergence cleaning method (e.g., [H] 02]) and a variant of the constrained 



transport method (e.g., [101 SSI IH3] ) , which are both discussed in detail in section 4.5 
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4.I. Evaluation of magnetic field expressions 

The MHD code stores the values of both the primitive B-field vector, 5* (Bvec in the 
code), and the evolved conservative field B^ (Bcons in the code) in the frame of the 
Eulerian observers. For analysis purposes we include options to compute the magnetic 
field in fiuid's rest frame 

,» = E^, (20) 

a 

b'=^ + WiB^v,)[v^-^-) , (21) 

1^2 



b' = ^ + iB\? . (22) 



(23) 



4-2. Reconstruction 



In a finite- volume scheme, one evaluates fiuxes at cell faces by solving Riemann problems 
involving potentially discontinuous hydrodynamic states on either side of the interface. 
To construct these Riemann problems, one must first obtain the fluid state on the left 
and right sides of the interface. The fluid state is known within cells in the form of 
cell-averages of the hydrodynamical variables. The reconstruction step interpolates the 
fluid state from cell averaged values to values at cell interfaces without introducing 
oscillations at shocks and other discontinuities. It is possible to reconstruct either 
primitive or conserved fluid variables, however using the former makes it much easier 
to guarantee physically vahd results for which the pressure is positive and the fluid 
velocities sub-luminal. 

By assuming an orthogonal set of coordinates, it is possible to reconstruct each 
coordinate direction independently. Thus, the fluid reconstruction reduces to a one- 
dimensional problem. 

We deflne t^j^i/2 to be the value of an element of our conservative variable state 
vector U on the left side of the face between Ui = U{xi, y, z) and f/j+i = [/(xj+i, |/, z), 
where Xi is the i^^ point in the x-direction, and t^/^i /2 the value on the right side of the 
same face (coincident in space, but with a potentially different value). These quantities 
are computed directly from the primitive values on the face. 

4-2.1. TVD reconstruction. For total variation diminishing (TVD) methods, we let 

where F{Ui) is a slope-limited gradient function, typically determined by the values 
of t/j+i — Ui and Ui — Ui-i, with a variety of different forms of the slope limiter 
available. In practice, all try to accomplish the same task of preserving monotonicity and 
removing the possibility of spuriously creating local extrema. GRHydro includes minmod, 
superbee [HI], and monotonized central [SS] limiters. TVD methods are second-order 
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accurate in regions of smooth, monotonic flows. At extrema and shocks, they reduce to 
first order. 

4-2.2. PPM reconstruction. The original piecewise parabolic method (oPPM, [75]) uses 
quadratic functions to represent cell-averages from which new states at cell interfaces are 
constructed. The fluid state is interpolated using a fourth-order polynomial. A number 
of subsequent limiter steps constrain parabolic profiles and preserve monotonicity so that 
no new extrema can form. The version implemented in GRHydro includes the steepening 
and flattening routines described in the original PPM papers, with a simplified flattening 
procedure that allows for fewer ghost points [221 ES]- The original PPM always reduces 
to first order near local extrema and shocks. The enhanced PPM (ePPM) maintains 
high-order at local extrema that are smooth [3S1 175] . 

4.2.3. ENO reconstruction. Essentially non-oscillatory (ENO) methods use a divided 
differences approach to achieve high-order accuracy via polynomial interpolation [7S1 [77] 
without reducing the order near extrema or shocks. In the third-order case, two 
interpolation polynomials with different stencil points are used. Based on the 
smoothness of the field, one of the interpolation polynomials is selected. 

4.2.4. WENO reconstruction Weighted essentially non-oscillatory (WENO) recon- 
struction [79] is an improved algorithm based on the ENO approach. The drawback 
of ENO methods is that the order of accuracy is not maximal for the available number 
of stencil points. Furthermore, since a single stencil is selected, ENO reconstruction is 
not continuous. In addition, the large number of required if statements make ENO meth- 
ods unnecessary slow. WENO reconstruction, on the other hand, attempts to overcome 
all these drawbacks. The essential idea is to combine all possible ENO reconstruction 
stencils by assigning a weight to each stencil. The weight is determined by the local 
smoothness of the flow. When all weights are non-zero, the full set of stencil points 
is used, and the maximum allowable order of accuracy for a given number of stencil 
points is achieved. When the flow becomes less smooth, some weights are suppressed, 
and a particular lower order interpolation stencil dominates the reconstruction. Using 
the same number of stencil points as third-order ENO, WENO is fifth-order accurate 
when the flow is smooth, and reduces to third order near shocks and discontinuities. 

We implement fifth-order WENO reconstruction based on an improved version 
presented in ^B]. We introduce three interpolation polynomials that approximate t/j^^ 
from cell-averages of a given quantity f/j: 

f^S/2= ^f^-2 - ^f/.-i + y t/. , (25) 

U-^l/2 = -\u^-l + \u^-\u.+l , (26) 

U-^y2= lu^ + \u.+,-\u.+2. (27) 



GRHydro: A new open source GRMHD code for the Einstein Toolkit 11 

Each of the three polynomials yields a third-order accurate approximation of U at the 
cell-interface. By introducing a convex linear combination of the three interpolation 
polynomials, 

f^^+i/2 = w'U^;\/2 + ^'U^;l/2 + ^'^m/2 , (28) 

where the weights w^ satisfy J2iW^ = 1; it is possible to obtain a fifth-order 
interpolation polynomial that spans all five stencil points {f/i-2, • • • , ^^1+2}- The weights 
w^ are computed using so-called smoothness indicators /3*. In the original WENO 
algorithm [79], they are given by 

13^ = -(4m2_2 - 19ui.2Ui-i + 25m2_^ + llui.2Ui - SlUi^iUi + IOm^) , (29) 
/32 = -(4m^_^^ - ISui^iUi + 13u^ + 5Mi_iMi+i - 13uiUi+i + 4m^+^) , (30) 
(3^ = -{lOuf - 3lUiUi+i + 25u^^i + llUiUi+2 - l9ui+iUi^2 + ^u'^+2) ■ (31) 



Note that, in this section only, /3* refers to the smoothness indicators (29) and not to 
the shift vector ([l]). The weights are then obtained via 

w' = —. ^ ^ , with w' = -. — '-—, , (32) 

where 7* = {1/16,5/8,5/16}, and e is a small constant to avoid division by zero. 
Unfortunately, the choice of e is scale dependent, and choosing a fixed number 
is inappropriate for cases with large variation in scales. The improvement made 



by [HS] overcomes this problem by modifying the smoothness indicators (29)-(31). 
Following [HS], we compute modified smoothness indicators via 

pi =pi + e\u^\+6 , (33) 

where \U'^\ is the sum of the Uj in the i-th stencil, and 6 is the smallest number that the 
chosen floating point variable type can hold. Now, the smoothness indicators depend 
on the scale of the reconstructed field. In our case, we set 

P' = (3' + ei\u^\ + l) , (34) 

with e = 10~2^ and use these values instead of /3* in (29)-(31). 



4-2.5. MPS reconstruction Monotonicity-preserving fifth-order (MP5) reconstruction 
is based on a geometric approach to maintain high-order at maxima that are smooth. In 
contrast to WENO or ENO, MP5 makes use of limiters to a high-order reconstruction 
polynomial similar to PPM to avoid oscillations at shocks and discontinuities. The key 
advantages of MP5 reconstruction are that it preserves monotonicity and accuracy, and 
is fast. In practice, it compares favourable to enhanced PPM and WENO in terms of 



accuracy (see Appendix B). MP5 reconstruction is carried out in two steps. In a first 



step, a fifth-order polynomial is used to interpolate cell- averages Ui on cell interfaces 



ttL . 



Uti/2 = (2t/i-2 - 13t/i-i + 47f/, + 27f/,+i - 3f/i+2)/60 . (35) 
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In a second step, the interpolated value U^^,2 is limited. Whether a limiter is applied 
is determined via the following condition. First we compute 

f/*^^ = Ui + minmod(f/,+i - f/„ a {Ui - f/i_i)) , (36) 

where a is a constant, which we set to 4.0, and where 

minmod(x,|/) = - (sign(x) + sign(?/)) min(|x|, \y\) . (37) 

A limiter is not applied when 

(t/ti/2 - U,){Ut,y2 - U'n <e\U\, (38) 

where e is a small constant, and \U\ is the L2 norm of Ui over the stencil points 
{Ui-2, ■ ■ ■ , Ui+2}- Note that the norm factor is not present in the original algorithm. We 
have added this to take into account the scale of the reconstructed field. Following [SD] , 
we set e = 10~^°, though in some cases, it may be necessary to set e = to avoid 
oscillations at strong shocks or contact discontinuities such as the surface of a neutron 
star. 



In case (38) does not hold, the following limiter algorithm is applied. First, we 
compute the second derivatives 

D; = U,^2 - 2Ui-i + f/, , (39) 

D^ = f/,_i - 2f/, + f/,+1 , (40) 

Dt = U,- 2f/,+i + f/,+2 . (41) 
Next, we compute 

Df' = minmod(4A° - D^, ADJ - D^ D^ Dj) , (42) 



Df t/2 = minmod(4A° - Dr ^ADr - D^ D^ D') , (43) 



'i+1/2 

-,M4 

where the four-argument minmod function is given by 

TmnTaod{w,x,y,z) = -{sign{w) + sign(x)) x (44) 

8 

|(sign(w) + sign(y))(sign(w) + sign(2;))| x 
min(|w|,|x|,||/|,|2;|)) . 
We then compute 

[/f^^ =U, + a{U, - t/,+i) , (45) 

U^"" =^(f/. + f/m), (46) 

U''^ = U^"" - lDt^/2 , (47) 

U'"" =U, + l{U.-U..-i) + lDf_{^,. (48) 

Using these expressions, we compute 

f/^in = max(min(f/i, f/,+i, f/*^^), mm{Ui, f/^^, f/^^)) , (49) 

f/^ax = min(max(f/„ f/,+i, f/^^^), max(f/i, f/^^, [/^^)) . (50) 

(51) 
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Finally, a new limited value for the cell interface f^/^1/2 is obtained via 

^L,limited ^ jjL^^^^ ^ minmod(f/mi„ - f/^^i/s, f/max " f^^i/2) • (52) 

To obtain the reconstructed value at the right interface U^^,2^ ^^^ values at 
Ui-2, ■ ■ ■ , Ui+2 are replaced by the values t/j+2, • • • , f^j-2- 

For a detailed description and derivation of the MPS reconstruction algorithm, we 
refer the reader to the original paper |8Uj . 

4-3. Riemann Solver 

We implement the Harten-Lax-van Leer-Einfeldt (HLLE) approximate solver [S21 EZ] ■ 
The more accurate Roe and Marquina solvers require determining the eigenvalues 
characterising the linearizing hydrodynamic evolution scheme, which is extremely 
resource intensive while not providing a decisive advantage in accuracy. In contrast, 
HLLE uses a two-wave approximation to compute the update terms across the 
discontinuity at the cell interface. With ^_ and ^+ the most negative and most positive 
wave speed eigenvalues present on either side of the interface (including magnetic field 
modes but not those associated with divergence cleaning subsystem, as discussed in 



section 4.5.1), the solution state vector U is assumed to take the form 

' U^ if < ^_ , 
U, if e_<0<e+, (53) 

5+ ! 



u 



U^ if 0>^^ 



with the intermediate state U* given by 

e+U^-^U^-F(U^) + F(U^) 
U* = T . (54) 

The numerical flux along the interface takes the form 

F(U) = ^+F(U^) - e-F(U^) + LU^^ - U^) ^ (55) 



t = min(0,^), 1+ = max(0,e+) • (56) 



where 



We use the flux terms (p5j) to evolve the hydrodynamic quantities within our Method 
of Lines scheme. 

Our calculations of the wave speeds is approximate, following the methods outlined 
in |13] to increase the speed of the calculation at the cost of increased diffusivity. The 
method overstates the true wavespeeds by no more than a factor of a/2, and then 
only for certain magnetic field and fluid velocity configurations. When computing the 
wavespeeds, we replace the full MHD dispersion relation by the approximate quadratic 
form (see 27 and 28 of jB]), 



2 7 2 



2 I 2 
■Va + Cs 



(57) 
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where we define the wave vector /c^ = {—uj,ki), the fluid sound speed Cg, the Alfven 
velocity, 



62 



62 



va 



(58) 
(59) 



I ph + b"^ V P^* 
the projected wave vector 

and the dispersion relation between frequency and (squared) wave number as 

Lo^ = kf.u'' = -uu^ + ku' , (60) 

kl=K^K^ = ul + g'"'k,k, . (61) 

The resulting quadratic may be written 



e' W\V' ~l)-V' -2i aW'v\V' - 1) + V'I3 



+ 



[aWef{V^ -l) + V^ (a'Y' - P'P' 



(62) 



where V"^ = v\ + c^(l — f^) and S, is the resulting wavespeed. Note that the indices are 
not to be summed over. Instead, we find different wavespeeds in different directions. 

When divergence cleaning is used to dissipate spurious numerical constraint 
violations that appear in the magnetic field, the characteristic wavespeeds include two 
additional {luminal) modes of the divergence cleaning subsystem. Since the divergence 
cleaning subsystem decouples from he remainder of the evolution system, its wavespeeds 
must be handled slightly differently 



which we discuss below in section 4.5.1 



4.4- Conservative to Primitive variable transformations 

In GRHD, converting the conservative variables back to the primitives is a relatively 



straightforward task, which can be accomplished by inverting (10) - (12) with all of 
the magnetic field quantities set to zero. GRHydro accomplishes the task through 
a ID Newton-Raphson scheme (summarised in section 5.5.4 of [22] and described 
in detail in the code documentation [Sni)- The scheme iterates by estimating the 
fluid pressure, determining the density, internal energy and the conservative quantities 
given this estimate, and using those in turn to calculate a new value for the pressure 
and its residual. The method works for any EOS, so long as one can calculate the 
thermodynamic derivatives dP/dp|e and dP/de|p. 

MHD adds several complications to the inversion, although the additional equation, 

E^ I ^. Instead, the primary difficulty is 



(13), is immediately invertible, yielding B^ 



that b^ cannot be immediately determined, and there is no simple analogue to the GRHD 
expressions that allow us to calculate the density easily once the pressure is specified. As 
a reminder, if we consider the (known) values of the undensitised conservative variables, 

— = pW , §' = — = ph*W%' - ab%' , (63) 

— = oh*W^ D* '-^^^"^ 

V7 



D 



T 



ph*W^ _ p* _ (aby 



D 
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the GRHD system (in which P* = P, h* = h, and 6^ = 0) allows us to define 

Q = f + D + P = phW^ , (64) 

and then determine the density as a function of pressure through the relation 



P= Q , (65) 

and thus the Lorentz factor and internal energy as well (see, e.g., ^OOj). In GRMHD, 
the most efficient approach for inverting the conservative variable set is often to use a 
multi-dimensional Newton-Raphson solver, with simplifications possible for barotropic 
EOS, for which the internal energy is assumed to be a function of the density only, 
eliminating the need to evolve the energy equation. Our methods follow very closely 
those of [03], particularly the 2D and ID\y solvers they discuss. 

We define a few auxiliary quantities for use in our numerical calculations. 
Unfortunately, it is impossible to construct a consistent notation that agrees with both 
the Einstein Toolkit release paper [22j and the paper that describes the conservative to 
primitive variable inversion scheme [53], so we choose to be consistent with the former 



here, noting the key differences in [Appendix A 



From the values of the conservative variables and the metric, we construct the 
momentum density 

S^ = -n,T\ = aTl , (66) 

whose spatial components are given by the relation Si = Si and its normal projection 
S given by 

S^ = {5'; + n^n'')S,. (67) 

Inside the Newton-Raphson scheme we make use of an auxiliary variable Q defined 



analogously to (64) by the expression 

Q = phW^ . (68) 

Since the inversion methods for general EOS and barotropic ones differ in some of the 
details, we discuss each in turn. 

4.4- i- General EOS The 2D Newton-Raphson approach implemented solves the 
following two equations for the unknown quantities Q and v"^, with all other terms 
known from the given conserved set: 

S' = S,S^ =v\B' + Qf-{S-Bf?^^, (69) 

S-n = -if + D)= -^(l + v')+^-^^Q-'-Q + P, (70) 

where all dot products are understood as four-dimensional when involving four-vectors 
and three-dimensional when not: e.g., 5 ■ B = S^B^^ which is equivalent to SiB^ since 
i?^ is a three- vector and thus S*^ = 0. 
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For a polytropic/Gamma-law EOS where P = [T — l)u and u = pe is the internal 
energy density, we may calculate P from the conserved variables and the current guess 
for Q and f ^ by noting that 

r-1 



P 



{1-v^)Q-DVi 



(71) 



The iteration updates Q and f ^ subject to the consistency conditions that < f ^ < 1 
and Q > 0, and a post-iteration check is performed to ensure that e > 0. 



For a more general EOS where (71) does not apply, and given the structure of the 



EOS interface in the Einstein Toolkit, it is easier to solve for the internal energy density 
u, and use it as a variable in a three-dimensional Newton scheme along with Q and v"^. 
First, noting that 

P = (1 - v^)Q - DVI - v^ - u , 



we may rewrite ( 70 ) as 



Sn 



B' 



:i+v') + 



{S ■ Bf 



2 ' ' 2 

In the Newton-Raphson steps, we first set p 



Q-^ -v'^Q-DVT 



v^ — u 



D\/\ — v"^ and solve 



u + p{p, u) = g(i - v^) - bVi^ 



(72) 



(73) 



(74) 



where the left-hand side will, in general, be a monotonic function of p and m, along 



with (69) and (73) 



Noting that the density depends only on the constant D and the value of v^ within 
the Newton-Raphson scheme, the partial derivatives of the pressure with respect to the 
Newton-Raphson variables are given by 



dP 

dP 

du 






dp 



(75) 



(76) 



We write these terms in this way, since the Einstein Toolkit EOS interface uses the 
pressure as a function of density and the specific internal energy e, P = P{p,e), and 
calculates partial derivatives against those variables, rather than p and u. 



4-4 ■ 2- Barotropic EOS For cases where the pressure and internal energy are functions 
of the rest mass density only, the EOS is barotropic (a polytropic P = Kp^ EOS is 
a special case of a barotropic EOS), we need a different inversion technique, since the 



quantity 5 • n in (70) requires knowledge of f, which is not evolved in these cases. 
Instead, the inversion uses (69) only, using it to eliminate the variable v"^ from the 



Newton-Raphson scheme by solving for v'^{Q): 



v\Q) 



Q\B^ + Qf 



(77) 
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In the special case of a polytropic EOS, we proceed by first solving for p{Q) through 
an independent Newton-Raphson loop over the equation 

,Q ^n^ii.'-^^). (78, 



Next, we use (77) and the fact that 



v^ = C-l, (79) 



to replace (69) by an expression given only in terms of Q and p{Q): 
Q = Q\B' + Q)W-Q\B^ + Qfv\ 

= Q^S^ + {S ■ Bf{B'' + 2Q) - l"^ - l"! Q\B^ + Qf . (80) 



When performing the Newton-Raphson step, all quantities in (80) are known a priori 



except the iteration variable Q and p((5), which depends upon it. The loop over Q 

(81) 



makes use of the derivative of (78), given by 
dp p 



dQ D^-iKp^-'^ - Q 

4-5. Divergence- free Constraint Treatment 

One of the main difficulties in numerical MHD simulations is the appearance of non- 



zero divergence of the magnetic field due to numerical errors, violations of (18) that 
would be interpreted as magnetic monopoles. A number of techniques have been 
designed to combat these, including "divergence cleaning" [1111121111], a method which 
introduces a new field that both damps and advects divergences off the grid, "constrained 
transport" algorithms [10] that balance out fiuxes exactly to maintain divergence-free 
magnetic fields to round-off accuracy, and magnetic vector potential methods that seek 
to achieve the same result [^ [22] • Constrained transport methods are often difficult 
to implement in simulations that employ mesh refinement, since maintaining balance at 
refinement boundaries is algorithmically complex. In the initial GRHydro MHD release 
we implement both divergence cleaning and constrained transport, both of which are 
detailed below. 

4-5.1. Divergence cleaning Divergence cleaning works by introducing a new field 
variable that both damps divergences and advects them off the grid through a hyperbolic 
equation modelled after the telegraph equation, driving numerical solutions towards 
zero divergence. Our implementation follows closely that of [H] |12], with some minor 
differences. The new field variable ip satisfies the evolution equation 

V,.(*F^V^^>)=«:n>, (82) 

a modification of Maxwell's equations ^ that reduces to the familiar form as ^ — > 0. 
Our parameter k, as we note below, determines the damping rate of the divergence 
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cleaning field ip, and incorporates the ratio of the parabolic to the hyperbolic damping 
speed dependence for the scheme, often denoted respectively by Cp and Ch in other works. 
The evolution equation for ip is given by the time component of (82). Noting 
that |73j 

1 



-iM^ 



W 



iu^'B'' - u^B 



''nt^\ 



(83) 



the first term on the left-hand-side of ( 82 ) yields, after some algebra and use of the fact 
that 5° 

1 



-aFOO = 0, 



V,.F^ 



and the second term 

Combining the two, we find 

dtip + di (aB' - i)p 



d^VlB' 



a^ 



-dtiP + P'd,tP 



v^(- 



-Ka 



d,j3']+^B^d, 



a 



v^y 



4) 



(85) 



(86) 



where we have grouped the derivative terms to respect the flux-conservative form. 

The spatial part of (82) is the modified evolution equation for the magnetic field. 
The two terms on the left hand side yield, respectively, 

{dty^B^ + d,^ [{av' - /3') B^ - (av^ - (3^) B'] } ,(87) 



V^F 



fj-j 






a^ 



(3Wti^+{a'Y^ -/3'(3ndi^ 



Using (86) to eliminate the dfip term, switching over to i3* as the primary magnetic field 
variable, and combining derivatives to produce a fiux-conservative form, we find after 
some algebra. 



dtB^ + d. 



[av 



(89) 



which reduces to the standard evolution equation for ■?/'—)• and diB^ — )■ 0. To evaluate 
the right-hand-side, we make use of the identity 

■1 



di{^l' 



-v^7''r^' 



kl 



^ 



'-i^'i'^'dau-i'S'dau 



(90) 



The characteristic velocities of the evolution system for MHD without and with 
divergence cleaning differ. In the former, the largest-magnitude wave speed is the 
fast magnetosonic wave speed. The inclusion of divergence cleaning introduces two 
additional modes, corresponding to eigenvalues for the evolution of the divergence 
cleaning field ip and the longitudinal component of the magnetic field, (i.e., the case 
where ^ = j inside the term in brackets in (90)), that in the fiat spacetime case decouple 

1. In our scheme, both modes 



from the remaining seven eigenvalues of the system 

have characteristic speed equal to the speed of light. Our system corresponds to setting 
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the hyperbolic damping speed Ch of the scheme to be the speed of hght in the notation 
of Newtonian divergence cleaning methods [IT]. This is in accordance with standard 
practice in both Newtonian and relativistic calculations to choose a characteristic speed 
for divergence cleaning that is as fast as allowed without violating causality, or in the 
case of Newtonian calculations, the Courant condition. We note that for the luminal- 



speed modes, the HLLE flux formula, (55) reduces to the local Lax-Friedrichs form when 
evaluated on a Minkowski background, 

F{U) = ^ [F(f/^) + F(f/^) - iiU"^ - f/^)] , (91) 

where ^ = max(^_|_, |,^_|). Test runs performed with the HLLE Riemann solver, 
but without the luminal-speed velocities often develop large spurious oscillations, 
particularly in cases where strong rarefactions are present, e.g., the cylindrical blast 
wave evolution described in section 15.31 and the rotor test described in section 15.41 

In our current implementation, we use the hyperbolic divergence cleaning technique 
only for flat spacetime tests. In this case the divergence cleaning subsystem decouples 
from the remainder of the MHD evolution system. In curved spacetime this decoupling 
is not as obvious and we opt to use constrained transport techniques instead, which we 
describe in the following. 

4-5.2. Constrained transport As an alternative to divergence cleaning, and to simplify 
comparison of the performance of the new code with existing GRMHD codes, we also 
implemented a variant of a constrained transport scheme [93] • In constrained transport 
schemes one carefully constructs a numerical update scheme such that the divergence free 



constraint ( 19 ) is conserved to numerical round-off accuracy. Rather than implementing 
the original scheme proposed in [nSHHS], which relies on a complicated staggering of the 
magnetic field components, we employ the simplified scheme called "fiux-CT" described 
in |ini SSI E3], a two-dimensional version of which is used in the HARM code [44j. The 
scheme uses cell-centred values of the magnetic field. For completeness of presentation 
we reproduce the basic description of the scheme found in |33], but refer the reader to 
the original literature for more details. 

In constrained transport schemes the induction equation ( [90) ) is written in terms of 
the electric field £ at the edges of each face of a simulation cell (see figure |l|. 

Employing Ampere's law, the time derivative of the face-averaged magnetic field 
component B^ is given by 

-^AyAz= -Sl^^^^^_^Ay-S^^y^^^,Az + Sl^^^^^^^Ay (92) 

+ ^1.-1,;^^^' 

with analogous equations for the other magnetic field components. In the ideal MHD 



approximation, the electric field components in (93) can be expressed in terms of the 
fluxes B'^v^ — B^v'^ for the magnetic field B given in ([l6J). Specifically, we use the 
numerical fluxes of the induction equation to calculate these electric field components. 
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Figure 1. Sketch of a simulation cell showing the location of face centred magnetic 
field and edge centred electric fields. 



Since the scheme of [13] evolves the cell-averaged magnetic field values, it computes 
the change in the cell-centred fields as the average change in the face centred fields 



?5m 

dt 



~dt 



+ 






(93) 



which is the actual evolution equation implemented in the code. 

In this formalism, the conserved divergence operator is given by 

I i+l A:+l 

(V ■ i3) .^1 ,,._^i ,,_^i = — — 2^ 2^ [Bi+i,j',k' - ^i,j',k' 



4Ax V •!,, I, 

J =j k'=k 



1 



i+1 fc+1 



4^ z2 z2 [^l,j+i,k' ^l,j,k' 

y i'=i k'=k 
1 «+l J'+l , ^ ^ . 

aT ^ ^ [^i',j',k+l " ^i',j',k) ' 



4A2^ , 



(94) 



which can be interpreted as a cell corner-centred definition of a divergence operator. 



5. Tests 



5.1. Monopole tests 

In order to test the divergence cleaning formalism and to explore its properties, we 
implement a number of tests that initialise the MHD configuration with numerical 
monopoles. Each case uses a uniform-density, uniform-pressure fiuid which is initially 
at rest, and we assume Minkowski spacetime. The ambient magnetic field is assumed to 
be zero. The initial magnetic configurations include "point" monopoles, where B^ ^ 
is set at a single point in the Center of the domain, and Gaussian monopoles for which 

e , r<nG , . . 

0; r>i?G, ^ ^ 
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where Rq is the radius of the compactly supported monopole. In the supplementary 
material we present additional results on the performance of the algorithm when dealing 
with high frequency constraint violations. 

In Figs. |2| we show the evolution of Gaussian and 3-dimensional alternating 
Gaussian monopole initial data. In each case, the grid is chosen to be 200^, spanning 
a coordinate range from —2 to 2 in each dimension, and we set Rq = 0.2. The 
simulations are performed using second-order Runge-Kutta (RK2) time integration 
and TVD-based reconstruction with an monotonized central limiter and a Courant- 
Friedrichs-Lewy (CFL) factor of 0.25 (i.e.. At/ Ax = 0.25 here). The fluid, assumed to 
follow a F = 5/3 ideal-gas law, i.e. to follow the condition P = (T — l)pe, was set to an 
initially stationary state with density and internal energy given by p = 1.0 and e = 0.1, 
respectively. We vary the divergence cleaning parameter k that appears in the driving 



term in (82), choosing values n = 1, 10, and 100. 

There is a markedly different evolution of the divergence of the magnetic field in 
time, following a predictable pattern. 

The choice k = 1 yields the slowest damping rate, which allows the wave-like 
behaviour of the divergence cleaning field to radiate away the divergence of the -B-field 
away from its original location as it damps downward. By contrast, k = 100 yields a 
much stiffer system of equations, resulting in a relatively slow damping rate, particularly 
for the lower-frequency terms in the initial data, with very low divergence in the magnetic 
field spreading through the numerical grid. The intermediate case, k = 10, yields the 
closest analogue to a critically damped system in that the amplitude of the divergence 
decreases most rapidly, while errors are efficiently radiated away across the grid. In 
general, we would recommend values of k ~ 10 to be used as a default. 

5.2. Planar MHD Shocktubes 

Historically, the simplest tests for an MHD scheme are shocktubes on a static fiat-space 
Minkowski metric. They are important to demonstrate the ability of the new code to 
capture a variety of MHD wave structures. Despite being simple setups they do provide 
a stringent tests for the algorithm. For our one-dimensional shocktube tests, we set up 
"Left" and "Right" MHD states on either side of a planar interface, with parameters 
drawn from the five cases considered in [HS] and summarised in Table [I| These have 
been widely used by several groups to establish the validity of numerical MHD codes. 
These parameter choices include generalisations of familiar cases long used to test non- 
relativistic MHD codes, e.g., case "Balsaral", which is a generalisation of the Brio-Wu 
shock tube problem ^7\ . 

To validate the new code, we evolve each of these shocks in each of the coordinate 
directions (i.e., "Left" states for domains of the grid satisfying x < 0, y < 0, and 
2; < and "Right" states for x > 0, y > 0, and z > 0, respectively), finding excellent 
agreement (limited only by numerical precision) in each direction, as expected. We 
also evolve shocks with mid-planes oriented obliquely to the coordinate planes, choosing 



GRHydro: A new open source GRMHD code for the Einstein Toolkit 



22 



CQ 



> 





t = 0.0 : 


jii t = 0.25: 


t = 0.5 : 








1 *^ : 


; . t = 0.75 : 


t=1.0 : 

_ • _ 

i 

> 


t = 1.25 : 


i 
i 




1 


; t = 1.5 : 

• 

! 


t = 1.75 : 


t = 2.0 : 


;■■■ K = 100: 

— 7C = 10^ 

; ... K = i - 


!■■ K = loo: 

— K = 10^ 
--■ K = \ - 


= ■- K = 100: 
: _ K = 10 ^ 

... K = \ 



-10 12 

X 

Figure 2. Behaviour of our divergence cleaning scheme, demonstrated by monopole 
damping and advection for a magnetic field with an initially Gaussian B'^ profile of 
radius Rg = 0.2. We show results for different values of the divergence cleaning 
damping parameter k. For k = 100 (dot-dashed blue curve), the system is quite stiff 
and damping is substantially slower than in other cases and more of the divergence 
remains localised near the initial position, while for n — 1 (dashed red curve) the 
damping is only slightly faster but spreads more rapidly away across the grid. For 
K = 10 (solid green curve), we find a nearly ideal choice of rapid damping and advection 
of the error away from the source, and recommend this value as a default for general 
calculations. 



x + y = ("2D diagonal") and x + y + z = ("3-d diagonal") as the mid-planes. 
Whereas shocks in any of the coordinate planes initially satisfy the divergence free 
constraint to roundoff error, and maintain this state indefinitely given the symmetry 
in the setup, diagonal shocks yield non-zero divergence of the magnetic field across the 
shock front due to finite resolution effects. The simplest way to understand this is to 
realise that for a 1-dimensional domain (say along the x direction), the x component 
of the magnetic field B^ is necessarily constant in space and time as a consequence of 
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the induction equation ( 90 ) . The slab symmetry of the test on the other hand imphes 
that all fields, in particular the B^ and B^ field, are independent of y and z. These 
two effects imply that diB^ = to roundoff precision initially and for all times, since 
the finite volume scheme preserves this property. Diagonal shocks, on the other hand, 
break this symmetry by having field components depend on two coordinates introducing 
divergence constraint violations of order of the truncation error. For non-flat geometries, 
the non-linear coupling of GRMHD equations would generate violations on the order 
of the truncation error in all directions. This multidimensional test allows us to verify 
that these constraint violations are controlled by employing divergence cleaning. 

In all cases, the coordinate origin was chosen so that no grid cell was centred directly 
on the shock front, so that each cell was unambiguously located in the left or right states. 
In each run we used RK2 time integration with a CFL factor of 0.8 for 1-d shocks and 
0.4 for 2D diagonal cases, with Ax = 1/1600 and Ax = Ay = 1/(800^2) for the 1-d 
and 2D diagonal cases to ensure equal resolution in the shock direction. For all results 
shown below, we use the HLLE Riemann solver and TVD reconstruction of the primitive 
variables with a monotonized central limiter. Boundary conditions are implemented by 
"copying" all data to grid edges from the nearest points in the interior located on planes 
parallel to the shock front. Divergence cleaning is turned on for the 2D case with k, = 10. 
Table [T] lists the parameters describing each shock test that is performed, and we note 
that all of them assume an ideal-gas law EOS with the given value of F, here F = 2 
for the "Balsaral" shock test and F = 5/3 for the others. These tests are also included 
as options in the code release within the GRHydro_InitData thorn. Exact solutions for 
each case are computed using the open-source available code of [HHl EH] • In all cases, 
we see very good agreement with the expected results that are comparable to results 
of similar codes in the literature [131 ESI EOl [Z31 ESI ESI ttSSl HDI]. We conclude that 
divergence cleaning does not interfere with generating physically meaningful GRMHD 
evolution results in the presence of shocks, and that the new code performs as well as 
other codes based on the same numerical techniques {89j . 

In figure |3] we show results for the relativistic generalisation of the Brio & Wu 
shock test develop in [HS]- The initial shock develops into a left-going fast rarefaction, 
a left-going compound wave, a contact discontinuity, a right-going slow shock and a 
right-going fast rarefaction. We find that the code captures all elementary waves and 
is in good agreement with the exact solution of |13] . In panel (g) of [3] we show the 
constraint violation as measure by a 2"^^^ order centred finite difference stencil for diB^. 
Figure 111 displays the results of the relativistic MHD collision problem (problem 4 of [nS] 
in which two very relativistic (Lorentz factor 22.37) streams collide with each other. We 
find deviations from the exact solution that are on the same level as in 

5.3. Cylindrical Shocks 

A more stringent multidimensional code test is provided by a cylindrical blast wave 
expanding outward in two dimensions. We take the parameters for this test problem 
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Figure 3. Evolution of the Balsaral shocktube, performed with divergence cleaning, 
showing ID and 2D diagonal cases as solid (blue) and dashed (red) curves, respectively. 
From top to bottom, we show the density p, gas pressure P, normal and tangential 
velocities v^ and v^ , the Lorentz factor W , the tangential magnetic field component 
B^ , and the numerical divergence of the magnetic field for the 2D case (for ID shocks, 
the divergence is uniformly zero by construction). Parameters for the initial shock 
configuration are given in Table [T] All results are presented at i = T^cU given by the 
final column of Table fll The results agree well with those reported in [JSl [HS] with 
most plots being indistinguishable and only the Lorentz factor plot showing a slight 
overshoot on the left-hand side of the shock. A more detailed description of the test 



setup and parameters can be found in the main text in section 5.2 
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Figure 4. Evolution of the Balsara4 shocktube case, performed with divergence 
cleaning, with all conventions as in figure |3] The results code reproduces the exact 
result very well, with deviations from the exact solution similar to those found inn |96j . 
A more detailed description of the test setup and parameters can be found in the main 
text in section [ 
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Table 1. Shock tube parameters. For each of these shock tube tests, originally 
compiled in [SB], we list for both the left "L" and right "i?" states the fluid density p, 
specific internal energy e, fluid 3-velocity v = u*, and tangential magnetic field B*, as 
well as the (uniform) normal magnetic field magnitude B", the adiabatic index T of 
the ideal-gas law equation of state, and the time Tj-^i at which results are plotted. Our 
values are quoted for shocks moving in the x-direction, i.e., the shock front is oriented 
in the y — z plane, so -B* is equivalent to {B'^ , B^) and S" = i?^. For diagonal shocks, 
all quantities are rotated appropriately along the diagonal oi x — y plane. Results for 
cases marked by an asterisk* are shown in the supplementary material only. 
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from |1U2] . The density profile is determined by two radial parameters, Tin and rout, and 
two reference states, such that 



p{r) 
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(96) 



Pout ; ?" ^ ''"out ! 

with an equivalent form for the pressure gradient. The initial fluid velocity is set to 
zero, and the initial magnetic is uniform in the domain. We use the following shock 
parameters: 



P. 



1-2 



10 



-4. 



Pin = 1.0, 



(97) 



-- 0.8, Tout = 1.0; Pin = 10"", Pout 

:3 X 10"^; B' = (0.1,0,0) . 

All tests use F-law equation of state with adiabatic index F = 4/3. We use a 200 x 200 x 8 
grid, spanning the coordinate range [—6, 6] in the x and y-directions, which is the setup 
used in |102j . To verify that our results are indeed convergent with resolution and 
to compare to [51], we run two more simulations using 400 and 500 points as well. All 
tests shown use the divergence cleaning scheme with k = 5, HLLE Riemann solver, TVD 
reconstruction using a monotonized central limiter. No explicit dissipation is added to 
the system of equations. 

This test is known to push the limits of many combinations of MHD evolution 
techniques, and often fails for specific combinations of reconstruction methods and 
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Rieinann solvers [OZl E2]. In our tests we only explore a small range of settings: we 
verify that the test succeeds when we swap TVD reconstruction with a 2^^^^ order ENO 
scheme or use constraint transport instead of divergence cleaning. In the later case we 
add Kreiss-Oliger dissipation |103j of order 3 and strength parameter e = 3 [lU3j to the 
magnetic field variables to stabilise the system. 

We find this test to be among the most sensitive one of our tests. Including the 



two light-like modes of the divergence cleaning field (see section 4.5.1) is crucial to 
obtain correct results. When these are not included, the code crashes quickly due 
to the growth of spurious oscillations in the magnetic field in the rarefaction region 
that develops behind the shock, particularly along the diagonals where the shock front 
expands obliquely to the grid (as was seen in tests of a several other codes [HSl 1^1^). 



We note that implementing the local Lax-Friedrichs flux (see (91 )) instead of HLLE also 
stabilises the evolution. 

Two-dimensional profiles are shown in figure |5] for P, VFLorentz, B^, and B^, along 
with numerically determined magnetic field lines, at t = 4. These results agree well 
with those shown in figure 10 of |102] . One-dimensional slices along the x— and y— axes 
for the rest mass density, gas pressure, magnetic pressure and Lorentz factor at t = 4 
are shown in figure [6] for three different numerical resolutions. Our results show good 
qualitative agreement with those presented previously by |57] and particularly by [55] . 
once one accounts for rescalings associated with different choices for the initial shock 
parameters. We find that the benefit of increasing the resolution is largest in the region 
within the shock front, which is consistent with results obtained by other groups [5711^ . 

5.4- Magnetic rotor 

As a second two-dimensional test problem we simulate the magnetic rotor test first 
described in [H5] for classical MHD and later generalised to relativistic MHD in [101] . 
The setup consists of a cylindrical column (the rotor) of radius r^^ = 0.1 and density 
Pin = 10 embedded in a medium of lower density pout = 1-0. Initially the pressure inside 
the column and in the medium are equal Pin = Pout = 1 and the cylinder rotates with 
uniform angular velocity oi Q = 9.95 along the cylinder axis, so that the fluid 3- velocity 
reaches a maximum value of Vmax = 0.995 at the outer edge of the cylinder. The fluid 
outside the cylinder is initially at rest. The equation of state used is a F-law, with 
r = 5/3. The cylinder is threaded by an initially uniform magnetic field of magnitude 
B^ = 1.0 along the x direction, covering the entire space; cylinder and exterior. We 
use TVD reconstruction with the minmod limiter, the HLLE Riemann solver, RK2 
time stepping with CFL factor 0.25. The magnetic field is evolved using the divergence 
cleaning technique using a damping factor k = 5.0. 

At the beginning of the simulation, a strong discontinuity is present at the edge of 
the cylinder since we do not apply any smoothing there. During the simulation, magnetic 
braking slows down the rotor while the magnetic field lines themselves are dragged from 
their initial horizontal orientation. At the end of the simulation, at t = 0.4, the field 
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Figure 5. Evolved state of the cylindrical explosion test problem, with parameters 
based on the test presented in |102] . At i = 4, we show the following quantities: 
top panels: Gas pressure P and Lorentz factor W, along with numerically determined 
magnetic field lines; bottom panels, B^ and B^ . The numerical resolution is Ax = 0.06. 
Profiles are almost indistinguishable from the solutions found in figure 10 of of ^102j . 
The very low pressure region outside of the shock from is slightly more extended 
in |102j . It is worth noting that |102j apply some numerical resistivity to control 
artifacts, while our divergence cleaning scheme has dissipative terms present in the 
divergence cleaning equations but do not add explicit resistivity, hence our results in 
low pressure regions might well differ slightly. 
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Figure 6. One-dimensional slices along the x— and y— axes for the evolved state of the 
cylindrical explosion test. Displayed are the rest mass density, gas pressure, magnetic 
pressure and Lorentz factor at t = 4. The red solid lines corresponds to a resolution of 
Ax = 0.06, the black dashed one corresponds to a resolution of Ax — 0.03 and the blue 
dotted on corresponds to a resolution of Ax = 0.024. These results may be compared 
to figure 7 of [57] and particularly figure 7 of [SS], noting that their parameters difi^er 
from ours, which were chosen to match the test as presented by |102j . particularly by 
the presence of a narrow transition region whose effects are visible as a secondary shock 
in the outer edges of the shock front. Comparing the plots, both codes agree on the 
overall structure of the result, with the two-dimensional pseudocolor plots in [5] being 
indistinguishable and the one-dimensional slices clearly showing the same structure of 
the shocks. A more detailed discussion of the test setup and results can be found in 
the text. 
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lines in the central region have rotated by nearly 90 degrees while at large radii the 
orientation of the field lines remains unchanged. 

Our results, shown in figure^ for a grid spacing Ax = 1/400, compare well with 
those presented in figure 5 of [FOlj and figure 8 of [69j. The density profile at the end 
of simulation is reproduced, showing slightly less noisy behaviour at the location of 
the expanding high density shell than the results reported in |101j . In addition, the 
magnetic field lines displayed in figure 5 of |101] are also reproduced very well. A slight 
over-density near the left and right corners of the (now low density) rotor is present in 
our simulations that is absent in |101j . 

Figure |8] depicts one-dimensional slices through the simulation domain at the end 
of the run. We show results for three different resolution, using 250, 400 and 500 
points to cover the domain. We find excellent agreement with previous work at similar 
resolutions 



5.5. Alfven Wave 

The propagation of low-amplitude, circularly-polarised Alfven waves has an exact 
solution that is useful for testing the performance of the MHD sector of the code [10] 
and provides a scenario for clean convergence tests, since the solution is smooth. The 
initial configuration is a uniform-density fiuid with a velocity profile given by 

v"^ = 0; t;^ = — f^y4ocos(/cx); f^ = —VAAQsm{kx) , (98) 

and corresponding magnetic field configuration 

B^ = const.; S^ = Ao5^cos(A;x); B' = AqB"" sm{kx) , (99) 

where we choose B^ = 1.0, an amplitude parameter Aq = 1.0, p = 1.0. A F-law EOS 
with r = 5/3 was used, as well as TVD as reconstruction method, a 2nd order Runge- 
Kutta (RK) [10411105] time integration with the Courant factor set to 0.2, and the HLLE 
Riemann solver, as well as divergence cleaning. 
The Alfven speed is given by the expression 
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(100) 



and the wavevector k = 2t\:/ L^ for one-dimensional cases. For two-dimensional cases, 
we rotate the coordinates so that wave fronts lie along diagonals of the grid. In all 
cases, periodic boundary conditions are assumed. Several values for the pressure have 
been used for testing purposes by other groups, spanning the range P = 0.1 [40j to 
P = 1.0 |106j . but here we choose P = 0.5, which yields the convenient result va = 0.5. 
In time, we expect the wave to propagate across the grid, such that if it travels an 
integer number of wavelengths we should, in the ideal case, reproduce our initial data 
exactly. In figure |9} we show the convergence results for our code for both 1-dimensional 
and 2-dimensional cases, finding the expected second-order convergence with numerical 
resolution when plotting the L2 norm of the difference between the evolved solution and 
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Figure 7. Evolved state of the magnetic rotor test problem, with parameters based on 
the test presented in llOlj . At i = 0.4, we show the following quantities: top panels: 
Density p and pressure P; bottom panels, magnetic pressure P„i = \ and Lorentz 
factor W , along with numerically determined magnetic field lines. Our results are in 
close agreement with those shown in |S71 [Ml llOlj . 



the initial data, evaluated along the x— axis (a full two-dimensional comparison produces 
the same result). We note that while the overall error magnitude appears larger than 
that shown for one-dimensional and two-dimensional waves in similar works, e.g., |1U6] . 
our results are shown after five periods, and we find comparable accuracy to |1U6] for 
the first cycle. 



5.6. Loop advection 

In this test, originally proposed in |lU7j and presented in a slightly modified form 
by [511 I1U6[ I1U8J . a region with a circular cross section is given a non-zero azimuthal 
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Figure 8. Evolved state of the magnetic rotor test problem, after t — 0.4. We display 
we display data for three different resolutions using 250, 400, and 500 cells each to cover 
the domain. The data was extracted along the x = 0.5 and y ~ 0.5 axis respectively. 
Our results are in very good agreement with figure 9 of [69J. 



magnetic field of constant magnitude (tlie "field loop"). The magnetic field is initialised 
to zero outside this region. The entire fluid configuration is given a uniform velocity 
and periodic boundary conditions are imposed. The analytical solution is for the field 
loop to move with constant velocity across the grid, with no magnetic field evolution. 
With Tc = y/x'^ + y"^ as the cylindrical radius we have an initial magnetic field 

o D I -^loopy/fci ^loopX l^ci ^ < -tiloop 1 /I ni \ 

""'■ "'^X 0; r > «,„, . 'l"^' 

In practice, this configuration results from the choice of the vector potential A = 
(0,0,m.8ix[0,Aioop{Rioop — ''^c)]), though we note that we here set the fi-field values 
explicitly in the initial data, rather than by finite differencing of the A-field, an option 
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Figure 9. Shown is the L2 norm of the difference between the evolved solution and the 
initial data of the Alfven wave test. Convergence is observed for the one-dimensional 
(ID, blue, dotted line) and two-dimensional (2D, red, dashed line) test, along with 
scaling expected for second-order convergence (solid line). In each case, as a function 
of N^, the number of points in the x-direction, we plot the L2 norm of the B-ficld 
along a line drawn in the x-direction through the centre of the computational domain, 
measuring the absolute difference between the initial and final values of B* after 5 full 
wave cycles for both the one-dimensional and two-dimensional cases. The resolutions 
used here were (16,32,64,128) and (20 x 15,40 x 30,80 x 60,160 x 120) points in ID 
and 2D respectively. 



also included in the code. 

Here, following |1U6] . we choose parameters Aioop 



10" 



P 



1.0, P = 3.0 for 



the initial magnetic field and hydrodynamic configuration. We consider two different 
models for the direction of the field loop: a "two-dimensional" case in which it is oriented 
vertically, and a "three-dimensional" case for which we rotate the direction of the tube 
so that it is oriented along the x — z face diagonal of our rectangular grid. For the initial 
velocity field of the two-dimensional configurations, we set v^ = 1/2, v^ = 1/24 and 
consider two different cases for the vertical velocity field parallel to the loop: one in which 
v^ = 0, so that the motion is perpendicular to the field loop, and one in which v^ = 1/24 
which additionally tests the ability of the conservative-to-primitive variable inversion 
scheme to maintain this quantity unchanged over time. For the three-dimensional case, 
we set v^ = 0.2^/2, v^ = 0.2, v^ = 0.1, and then rotate the tube and all vectors by 
7r/4 radians in the x — z-plane so that the loop runs along the x — z face diagonal. 
In all cases, velocity components are chosen to yield integer intervals over which the 
flux tube should propagate across the grid, T = 24 for the two-dimensional cases and 
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T = 5 for the three-dimensional case. Runs were performed using RK2 timestepping, 
TVD reconstruction with an monotonized central limiter, the HLLE Riemann solver, 
divergence cleaning active with k, = 10, and CFL factors of 0.4 for the 2-d cases and 
0.25 for the 3-d case. 



In the top panels of figure 10, we show the initial state of the flux tube of the 
two-dimensional case, with B^ on the left and the magnetic pressure Pm = &^/2 on the 
right. In the centre panels, we show the results of a low-resolution simulation, with 
Ax = 1/128, and in the bottom panels a high-resolution simulation with Ax = 1/256. 
Both simulations are evolved for one grid crossing time, T = 24, with v^ = 1/24. We find 
no substantial differences in the results after one period between v^ = and v^ = 1/24. 
While the observable decrease in the central magnetic pressure is the inevitable result of 
the simulation, we find the effect is significantly decreased by the increase in resolution. 
We find similar results for the three-dimensional loop, whose evolution is depicted by 



figure 11, and demonstrates comparable uniformity in the pressure and field strength 
for a given resolution. From top to bottom, we show the initial state and the states 
after one and two periods (T = 5 and T = 10, respectively). The initial pressure, which 
should remain invariant, is maintained throughout the body of the loop, but decreases 
in the central region and near the edges, as expected (this result can be compared with 
figure 4 of |106j . with which we find excellent agreement). 

5. 7. Bondi Accretion 

Bondi accretion is a steady-state solution for spherically-symmetric, adiabatic accretion 
onto a black hole. The Einstein Toolkit implementation of the general relativistic 
analogue of Bondi's original solution jl09j follows the methodology laid out in |110j . 
assuming a fixed background spacetime of a Schwarzschild black hole of mass Mbh- 
The solution is determined by the continuity equation 

M = ATTpur^ , (102) 

and the integrated Euler equation for a polytropic equation of state, 

Q = h\l-^-^ + e), (103) 

where " denotes a variable in Schwarzschild coordinates, h is the specific enthalpy, 
u = —Uir\ and Cs the speed of sound. Adding a radial magnetic field to the system 
does not alter the Bondi solution. Following [14j, we set the radial magnetic field by 
prescribing a primitive magnetic field vector of 

B- = ^^ , (104) 



where r = y/x'^ + y'^ + z^ is the coordinate radius and Bq controls the magnetic field 
strength. A magnetic field vector of this form by construction satisfies the divergence 



free constraint (19). 



A particular solution is chosen by specifying the BH mass Mbh, the mass accretion 
rate M, the location of the sonic point in Schwarzschild coordinates (i.e., the areal 
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Figure 10. Results of the two-dimensional advected loop test with u^ — 1/24. 
In the top panels, we show the initial state of B^ in greyscale, along with the 
(uniform) magnetic pressure. In the centre and bottom rows, we show the result 
of low (Ax = 1/128) and high (Aa; = 1/256) resolution simulations after one period 
has passed (T ~ 24) . The loss of magnetic pressure in the centre of the loop is greatly 
reduced in the higher-resolution simulation. 
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Figure 11. Results of the three-dimensional advected loop test with v^- ^ 0. In the 
top panels, we show the 2d z = slices through the initial state of B^ in greyscale, 
along with the (uniform) magnetic pressure. In the centre and bottom rows, we show 
the results after one (T — 5) and two (T = 10) periods, respectively. The flux tube 
formed by the moving flux loop is oriented diagonally in the domain, hence the slice 
is elliptical rather than circular. The grid spacing used is Ax — 1/128. Our result is 
extremely similar in how well it preserves the magnetic flux and pressure of the loop 
to the results shown in figure 4 of [106 . 
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radius) fg, the magnetic field strength Bq and the adiabatic index F of the equation of 
state. Table |2] hsts the values of the parameters used in our simulations. We explore 
various values for Bq and M but keep the sonic radius f = 8 Mbh and adiabatic index 
r = 4/3 fixed. 

By requiring the solution to remain smooth at the sonic point, the sound speed 
and thus velocity of the fluid at that point must satisfy Ug = Cs = MBH/(2fs). From 
this equation we find the polytropic constant K for the EOS and subsequently from the 
continuity equation the conserved quantity Q. The solution in the remaining domain 
is then determined by integrating outwards (or inwards) from the sonic point before 
converting to the desired coordinate system of the evolution. We choose to work in 
ingoing Kerr-Schild coordinates which are horizon penetrating and therefore allow us to 
model the flow of matter across the event horizon. 

To avoid numerical issues due to the singularity inside of the event horizon we use 
a transition function 



/ / 7rr 71 
tanh tan 



X(r) = <^ 2 L V VMbh 2 




r < Mbh , , 

(105) 

otherwise 



to damp magnetic field B^ — > xi?* and velocity v^ — )■ xf * to zero near the singularity. 
Similarly we reduce the mass density, pressure and internal energy density to values 
corresponding to our atmosphere level Patmo which is 10^*^ of the density at the event 
horizon. 

Our numerical setup follows [13], i.e. we simulate a rectangular box of size 
[0,11] X [0,11] X [—11,11] using resolution of 0.11 for the low resolution and 0.073 
for the high resolution run. Since the system is spherical, we used a 90 degree rotational 
symmetry around the z axis to reduce the computational costs. We use a HLLE Riemann 
solver, TVD reconstruction using a minmod limiter. We integrate in time using a 
third-order Runge-Kutta (RK3) integrator with Courant factor 0.25. For this test we 
employed our constrained transport scheme which we found to be more stable in the 
strongly magnetically dominated regimes explored by us in this test. After each RK 
substep we re-set the grid inside of f < 1 < Mbh and at the outer boundary to the 
analytic initial solution, providing Dirichlet type boundary data. 

We note that due to our choices of numerical timestepping method, stencil size 
and resolution grid points outside of the event horizon f = 2 Mbh are not affected by 
grid points inside of f = 1 Mbh- Thus our method is effectively excising the interior of 
r < 1 Mbh from the simulated region. 

To verify the accuracy of our code, we evolve the analytical initial data for 
T = 100 Mbh and compare the data along the z = 0.055, y = line to the analytic 



initial data. We display the result of this comparison in figure 12 We observe 2°^^ 
order convergence with increasing resolution in the region outside of the event horizon 
at r = 2 Mbh- Our results compare very well with those of [33], figure 6. 

To adequately test the stability of the steady-state solution, we repeat this test 
and vary both the accretion rate and the magnetic field strength. We also include a 
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Table 2. List of tests performed for the magnetised Bondi infall. The first columns 
shows the run name, the second the magnetic field strength parameter, the third the 
(unitless) mass accretion rate, the fourth the ration of magnetic pressure to gas density 
at the horizon, the fifth the resolution and the last column lists the maximum relative 
error in p outside of r = 2M for each test run. 



Run 


Bo 


M 


bVp 


Ax 


|(p/Poxact)- l|oo 


A 


0.0 


12.57 


0.0 


0.11 


9.70 • 10-^ 


B 


0.0 


12.57 


0.0 


0.073 


4.60 


10-^ 


C 


11.39 


12.57 


25.0 


0.11 


7.42 


10-2 


D 


11.39 


12.57 


25.0 


0.073 


3.92 


10-2 


E 


11.39 


25.13 


12.5 


0.11 


3.88 


10-2 


F 


5.7 


12.57 


6.26 


0.11 


2.21 


10-2 


G 


5.7 


25.13 


3.13 


0.11 


1.37 


10-2 



pure hydrodynainic simulation for each accretion rate. In Table [2| we list the set of 
parameters and their values that we vary for this test along the supremum norm of the 
relative error on the grid outside of the event horizon after t = 100 Mbh of evolution. 

We found that higher order reconstruction method such as PPM and methods 
employing divergence cleaning performed less well in this test. In the case of divergence 
cleaning we suspect that our inability to specify initial data with vanishing numerical 
divergence is partially the source of these difficulties. Even in the constrained transport 
case, we observe a certain amount of initial drift away from the analytic solution in the 
field values. Eventually the solution settles down and we the amount drift to decrease 
with increasing resolution. In the case of divergence cleaning, our imperfect boundary 
data serves as a constant source for constraint violations, possibly spoiling the result. 
Similarly the interpolation inherent in adaptive mesh refinement techniques seems to 
serve as a source of constraint violation and prevent us from simulating this system 
using AMR and either divergence cleaning or constraint transport. 

5.8. Magnetised TOV Star 

In this section we consider the stability of a stationary, poloidally magnetised spherically 
symmetric neutron star model, obtained via the solution of the Tolman-Oppenheimer- 
Volkhoff equation |lllj . We follow [IT] and set up a poloidal magnetic field on top of 
the independently solved fiuid configuration. This setup provides a more challenging 
test for the code than the previous ones since it probes both the GR and the MHD 
evolution. While the initial configuration is stationary, we evolve both the spacetime 
and the GRMHD equations to test the ability of the new code to maintain such a 
solution. The hydrodynamic initial data is a TOV star generated using a polytropic 
EOS with polytropic constant K = 100, adiabatic index F = 2.0 and initial central 
density of 1.28 ■ IO-^Mq. The fiuid is evolved with a T-law EOS and T = 2. As 
in [121, we introduce a poloidal field that we add on top of the TOV initial data. 
We choose the vector potential A^ = Ai,zu'^{l — po/p™'^^)"''' max (P — Pcut,0) where 
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Figure 12. Bondi inflow after IOOMbh of evolution compared with the semi-analytic 
solution for varying numerical resolution. We observe 2"^^ order convergence in the 
region outside of the event horizon at r = 2 Mbh ■ The red dashed line shows the 
errors in p compared to the analytical solution for the low resolution run rO with 
resolution 0.11 while the black solid line shows the error in the high resolution rim 
rl using a resolution of 0.073 scaled so that for second order convergence the curves 
coincide with each other. This inset shows a zoom in of the region around the event 
horizon where matter is moving fastest and the magnetic field is strongest, i.e. where 
we expect the largest errors. Our errors are almost identical to those found by |43j . 



zu = \/ [X — x^)2 + y^. Here Af, determines the strength of the initial magnetic field, Up 
the offset of the maximum of the magnetic field strength in respect to the maximum 
density pmax of the fiuid, and Pcut the pressure outside of which the magnetic field is set 
to zero. In this test we confine the magnetic field to the interior of the star by setting 
-Pcut = lO^^Patmo, much higher than the atmospheric pressure. The parameter Ah is 
set to unity, which, together with Up = 0, provides a magnetic field with a maximum 
strength of 8.5 x 10~^ in geometrised units which corresponds to 2.5 x 10^^ G. 

We place the TOV star with radius S.ISMq and gravitational mass I.4OM0 at the 
origin of a grid with 5 levels of mesh refinement with cubical extents 32OM0, I2OM0, 
6OM0, 3OM0, and 15M0. The extent of the finest level is chosen in order to cover the 
entire star. We perform three simulations with coarsest resolutions of 8Mq (rO). 4:Mq 
(rl), and 2M0 (r2). The corresponding resolutions of the finest grids are I.OM0 (rO), 
O.5M0 (rl) and O.25M0 (r2). We simulate for 6.9 ms using oPPM reconstruction, the 
HLLE Riemann solver, and evolve the magnetic field with constrained transport. We 
use 4th-order Runge-Kutta time integration with a Courant factor of 0.25. We set the 
F-Driver gauge condition [SH] parameter to 77 = 0.5 and add Kreiss-Oliger dissipation 
to the spacetime variables with a dissipation coefficient of e = 0.1. 



Figure 13 shows the fractional change of the maximum rest-mass density p for 
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Figure 13. Results of the magnetised TOV simulation. First panel: Normalised 
maximum rest-mass density Pmax/pmax.o- Second panel: L2-norm of the Hamiltonian 
constraint ||-ff||2 for the three resolutions rO, rl, and r2. Third panel: Absolute 
difference e in normalised maximum rest-mass density between low (rO) and medium 
(rl), and medium (rl) and high (r2) resolutions are shown (centre panel). Fourth 
panel: L2-norm of the normalised divergence of the magnetic field \ViB'^/{\B\ ■ 
Aa;~^)|2, where Ace is the finest grid spacing, for each of the the three resolutions rO, rl, 
and r2 and |-B|oo is the maximum magnetic field strength on the grid. The differences 
and the _L2-norm of the Hamiltonian constraint and the divergence of magnetic field 
are scaled for first order convergence. 
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all three resolutions, the differences between low (rO) and medium (rl), and medium 
(rl) and high (r2) resolutions in maximum rest mass density and the L2-norm of the 
Hamiltonian constraint | |-ff 1 12- The differences in rest-mass density and the Hamiltonian 
constraint are scaled for first order convergence. While the TOV solution is stationary, 
numerical perturbations excite oscillations in the central density of the star, commonly 
observed in simulations of TOV stars, which must converge away with increasing 
resolution. We can analyse the convergence properties of our simulations more carefully 



by considering the second panel of figure 13, which shows the absolute differences 



between low (rO) and medium (rl), and medium (rl) and high (r2) resolutions in the 
normalised maximum rest-mass density. The differences are scaled by a factor of 2 for 
first-order convergence. We observe between first and second order convergence, which 
is the expected behaviour, since the oPPM scheme reduces to first-order accuracy at 



the surface of the star. The third panel of figure 13 shows the L2-norm of the evolution 
of the Hamiltonian constraint ||-ff||2 for the resolutions rO, rl and r2. The medium 
and high resolutions are rescaled for first-order convergence. We again observe between 
first-order and second-order convergence, consistent with what is expected. The fourth 



panel of figure 13 shows the L2-norm of the normalised divergence of the magnetic field 
as a function of time for the resolutions rO, rl and r2. Again, we observe between first- 
order and second-order convergence. Note that the lines are constant, since constrained 
transport preserves the initial constraint violation to round-off. 

5.9. Rotating Core Collapse 

As a final test, we present results for the collapse of a rotating magnetised stellar core 
with simplified microphysics. We analyse a number of quantities to demonstrate the 
ability of our code to simulate magnetised stellar collapse in full GR. For this, we 
consider model Al of Shibata et al. |T2], who carried out simulations in axisymmetry. 
While our code is 3D, we can limit non-axisymmetric dynamics to even i and m (in 
terms of spherical harmonics) that are multiples of 4 by evolving only in one octant 
of the 3D domain and enforcing rotational symmetry on the two inner faces of the 
octant and reflection symmetry in the z-direction. We model the precollapse star as a 
r = 4/3-polytrope with a central density of pc = 10^*^ g cm~^ and a polytropic constant of 
K = 4.897- 10^''[cfy's]. This corresponds to an ultra-relativistic degenerate electron gas at 
electron fraction Ye = 0.5. We choose the initial rotation profile as u^u^ = ip"^ {flc — ^), 
where flc is the central angular velocity, Q is the angular velocity, and ipd is a constant 
describing the degree of differential rotation. Model Al of [19] is uniformly rotating 
ifd -^ oo) and has a ratio of polar to equatorial radius Rp/Re = 0.667. The initial 
data were generated using the RNS code |112j . During evolution, we employ a hybrid 
EOS |113t I114J . consisting of a two-piece polytropic cold pressure component Pp and a 
thermal pressure component Pth, 

P = Pp + Pth . (106) 
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Figure 14. First panel: Convergence of tlie gravitational wave signal A2 computed 
via the quadrupole formula as in J19j . and gravitational wave signal from [19^ for the 
corresponding model. Second panel: Central rest-mass density pc- Third panel: L2- 
norm of the Hamiltonian constraint ||^||2- Fourth panel: L2-norm of the normalised 
divergence of the magnetic field iWiB"^ /{\B\ •Aa;~^)|2. All four quantities are shown 
for the three resolutions rO, rl, and r2. The L2-norm of the Hamiltonian constraint 
and the normalised divergence of the magnetic field in panels three and four are scaled 
for first order convergence. The constraints initially exhibit second-order convergence, 
which drops to first order after core bounce generates a shock. The divergence of 
the magnetic field is kept constant by the constrained transport scheme, unless a 
change in the grid structure occurs. This causes the jumps in the time evolution of 
|Vi-B'/(|i?| • A.T-^)|2 since the newly created grids are populated using non-constraint- 
preserving operators. 
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Figure 15. Convergence of the maximum value of the magnetic field components B^ 
for the three different resolutions rO, rl and r2 in the rotating core collapse simulation. 
Differences between the three resolutions are most obvious for i?^ (top panel) , but also 
present for B^ (middle panel) and B^ (bottom panel). This is however to be expected, 
as differences in the maximum of the magnetic field components may arise from local 
turbulent behaviour and depend sensitively on resolution. We observe a similar overall 
evolution of the magnetic field maxima as in [TS]. We note that the differences in 
evolution for i3f,„^ and B^^^ result from the fact that we choose to plot the maximum 
value and not the maximum of the absolute value. As our domain is symmetric under 
rotations by 90 degrees about the z-axis, there is a sign flip between B'^ and B^ in any 
given quadrant of our domain, i.e. S^^x = ~-^min- 



Following ^19j, in the cold part, the adiabatic index jumps from Fi = 1.3 to F2 = 2.5 
at p > 2 X lO^^gcm"^ to mimic the stiffening of the nuclear EOS at core bounce. The 
thermal component of the hybrid EOS is a F-law with Fth = Fi. The specific internal 
energy is given by e = ep + Cth and consists of a polytropic and a thermal part. We 
induce collapse by following the prescription outlined in [12], that is we set the internal 
energy according to e = 3Kp^^^ and then determine the pressure as P = 3i^(Fi — l)p^/^. 
This reduces the pressure from its initial value by a factor of (1 — 3(Fi — 1)) (i.e. 10% 
for Fi = 1.3). Following [TU], we assume a poloidal initial magnetic field computed from 
the vector potential 



Aj, = AhW max 



P Pcut : 
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Table 3. Initial parameters and properties of model Al. 

Central rest-mass density [gcm-^] 

ADM mass [Mq] 

Baryonic mass [Mq] 

Ratio of rotational kinetic to potential energy 

Equatorial radius [km] 

non-dimensional angular momentum parameter 

Ratio of polar to equatorial radius 

Differential rotation parameter 

Angular velocity at the rotation axis 

Maximum initial magnetic field strength [G] 

Maximum initial ratio of magnetic and fluid pressure 

Ratio of magnetic to rotational kinetic energy 



with 77.6 = 1, Pcut = 10~^Pc and w = y/x'^ + y^ . The parameter A\, is chosen so that the 
maximum initial magnetic field strength is 7.2 x lO^^G. For convenience, we summarise 
the model parameters in Table [3j Our numerical setup initially consists of 4 levels of 
nested grids each refined by a factor of 2 in resolution. As the star collapses, and the 
central density increases, we progressively add 5 additional refinement levels, following 
the approach taken in |115j . The density thresholds are chosen to be [8. Ox 10^°, 3.2x 10^^, 
1.28 X 10^^, 5.12 X 10^^, 2.05 x 10^^] g cm"^. The final grid structure consists of 9 levels of 
nested grids with cubical extents [5669.8 km, 3061.7 km, 2438 km, 1559.2 km, 283.5 km, 
212.6 km, 144.7km, 59.1km, 23.6 km]. The extent of the finest level is chosen such that 
it covers the entire compact remnant (the "proto-NS") after core bounce. The inner- 
most refinement level for our fiducial setup (rl) has a linear resolution of 369.1m. We 
consider two more resolutions, rO with fine grid resolution 461.4 m and r2 with fine grid 
resolution 295.3 m. This allows us to study the resolution dependence of our results. We 
use vertex-centred AMR, constrained transport, and TVD reconstruction with the van 
Leer monotonized central limiter. While GRHyro includes higher-order reconstruction 
schemes, we use TVD in these simulations to match the numerical methods used in [12] 
as closely as possible. The time integration is performed using RK4 with a Courant 
factor of 0.2. The damping coefficient of the F-driver gauge condition for the spacetime 
evolution is set to r^ = 1/2. We add 5th-order Kreiss-Oliger dissipation to the spacetime 
variables with a dissipation coefficient e = 0.1. The atmosphere is chosen to be a factor 
10^° smaller than the central density. Shibata et al. [12] extract the gravitational wave 
signal from rotating core collapse and bounce via a quadrupole formula and define a 
wave amplitude 

A2 = h. - hz , (108) 

with the definitions of Ixx and Ixx given in [19] . Since the gravitational wave signal is 
a good tracer of the overall dynamics of rotating core collapse, we extract A2 from our 
simulations for comparison with |19j . 
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In figure [14], we show convergence for tfie central density pc, tlie gravitational 
wave signal A2, the L2-norni of the Hamiltonian constraint ||-ff||2, and the normalised 
divergence of the magnetic field \ViB^/{\B\ ■ Ax~^)|2. The Hamiltonian constraint in 



the third panel of figure 14 and the normalised divergence of the magnetic field in the 
fourth panel are scaled for first-order convergence. Before core bounce, second order 
convergence is achieved, but reduced to first order after bounce. Our numerical scheme 
reduces to first-order accuracy at shocks, contact discontinuities and local maxima, i.e. 
at the centre of the proto-NS, which dominates the Hamiltonian constraint violation. 
The gravitational wave signal A2, the central rest-mass density pc, the L2-norm of the 
Hamiltonian constraint and the normalised divergence of the magnetic field are shown 
for the three resolutions rO, rl, and r2. In the top two panels, we also show the 
gravitational wave signal and central density evolution as extracted and digitised from 
figures 2 and 17 of [IS]. We observe excellent agreement for the gravitational wave 
amplitude A2 with the results of pL9j. The central density evolution of [12] reaches a 
peak value of 6.5 ■ lO^^gcm"^ at bounce while we find a maximum bounce value of 
6.4 ■ lO^^gcm"^. The post-bounce evolution differs by an approximately constant offset 
of ^ 6%, however the features in the time evolution itself are very similar. We attribute 
the observed deviations in maximum rest-mass density evolution to differences in our 
numerical setup. 

We find between second-order and first-order convergence for the gravitational 
wave signal A2, the central density Pc, and the L2-norm of the Hamiltonian constraint 
||-ff||2, which is within the expected behaviour of the code. The convergence in 
both central density and maximum rest-mass evolution are less clean than for the 



Hamiltonian constraint due to their oscillatory nature. The fourth panel of figure 14 
shows the evolution of the L2-norm of the normalised divergence of the magnetic field 
|Vi-BV(|-B| • Ax^^)|2. The constrained transport scheme used in this test maintains 
the divergence present in the initial data to round-off precision. Constraint violations 
appear only when additional grids are added and filled by non-constraint-preserving 
interpolation. 



Figure 15 shows the maximum of the magnetic field components |-B*|max as a 
function of time for the three resolutions rO, rl, and r2. We observe resolution- 
dependent differences in the magnetic field evolution in the postbounce phase. This is 
particularly obvious for B^, but also visible in B^. This behaviour is expected, since the 
evolution of the magnetic field is very sensitive to resolution due to effects of turbulence 
and instabilities, most importantly the magnetorotational instability (MRI) [116^ 1117] . 
A similar behaviour for the resolution dependence of B^ is shown in figure 6 of [12]. Our 
results for the postbounce magnetic-field evolution differ quantitatively from [TU] , which 
is expected due to our different numerical setup. While our setup consists of different 
layers of mesh refinement with varying resolution throughout our 3D domain, [19j employ 
a uniform axisymmetric grid with a coordinate singularity at the rotation axis. 
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6. Conclusions 

We have presented the implementation and a set of rigorous tests of ideal general- 
relativistic magnetohydrodynamics (GRMHD) in the GRHydro code, which is part of 
the Einstein Toolkit. All tests and results presented in this paper can be reproduced 



with the code and input files available at http : //einsteintoolkit . org 



GRHydro is the first open-source code that provides for fully self-consistent GRMHD 
evolutions on dynamical general-relativist ic spacetimes. Its development includes 
many contributions from a wide array of collaborators involved with other relativistic 
codes, particularly the original version of the general-relativistic hydrodynamics code 
Whisky [32j, from which GRHydro started, and the publicly available fixed-background 
HARM code |11|, from which we have adapted a number of routines and techniques. 
Combined with other functionality from the Einstein Toolkit, the code may be used 
to perform GRMHD evolutions of a broad range of relativistic astrophysical systems. 
We expect these to primarily consist of magnetised isolated and binary neutron stars 
and collapsing stellar cores, but should the proper initial data be constructed, accretion 
disks and tori, white dwarfs, and a number of other configurations could be simulated 
as well. 

Our implementation of the equations of GRMHD also features the ability to 
eliminate spurious divergence in the magnetic field through two different methods: 
hyperbolic divergence cleaning and constrained transport. While these are clearly 
sufficient for many problems, one of the responsibilities in implementing a community- 
based code is to incorporate techniques developed by other groups to solve similar 
problems. Along these lines, we foresee work in the future to incorporate vector 
potential treatments in which the 5-field is computed from the values of an underlying 
vector potential A that serves as a fundamental evolved variable. Such a scheme has 
recently been implemented by Etienne et al. [HSl 02] and Giacomazzo et al. |118] . The 
former work also included appropriate electromagnetic gauge conditions for evolutions of 
compact object binaries involving neutron stars and other relativistic systems of interest. 

The GRHydro code now includes several improved higher-order techniques for 
reconstruction of states at cell interfaces for use in high-resolution shock capturing, 
including three that were described in the original Einstein Toolkit code paper [22], 
TVD, PPM, and ENO, as well as three new methods that have been implemented and 



released since, ePPM [Ml EH], WEN05 [THl EH] and MPS [80j (see|AppendjxB|. 

In the near future, we will extend the new GRMHD version of GRHydro to 
the generalised multipatch infrastructure of [36], which will be released as part of 
the Einstein Toolkit and will enable simulations on a combination of Cartesian and 
curvilinear grids. Further work will be directed towards implementing improved 
approximate Riemann solvers for GRMHD. While the Roe |119] and Marquina 1120(1121] 
solvers implemented in GRHydro for purely hydrodynamic evolutions are likely to 
be too computationally expensive, requiring an eigenmode breakdown of an 8 x 8 
matrix system at every grid cell, the HLLC solver, which recovers the contact wave 
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dropped by the HLLE approximation, may be an attractive and computationally 
efficient improvement [7H 1122] . Additionally, we will introduce improved techniques for 
performing conservative-to-primitive variable conversions to handle cases in which the 
current routines are known to fail, particularly those where the Lorentz factor becomes 
very large or the fluid becomes dominated by magnetic pressure |123] . 

In the long term, we expect the GRHydro code and the Einstein Toolkit to 
evolve into a toolkit not just for numerical relativity, but for the broad computational 
relativistic astrophysics community. Towards this end, we will include additional physics 
capabilities in future versions of the Toolkit. These will include better treatments of 
finite-temperature and composition dependence, nuclear reaction networks, radiation 
transport of neutrinos and photons, and improvements over the current ideal GRMHD 
approximation, following the most recent developments in resistive GRMHD |124] 1125] . 
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Table Al. Overview of differences in variable notation between HARM and Valencia. 



Quantity 


Valencia/ET 




HARM 




Pressure 


P 




P 




Energy density 


ph = p{l + e) + P 


w 




3- velocity for a normal observer 


, u' 13' 

7 a 
W 




v' 




Lorentz factor 




7 




Projected velocity 


Wv' 




u' ^ jX - 




Momentum density for a normal 


S, ^ aTf, 




Qm 




observer 










Projected momentum density 


Sfj, = (5'^ + n^n 


")S. 


Qm 




Auxiliary Con2prim variable 


Q = phW^ 




w 





visualisation package [nJlll26j . All other figures were generated with the Python-based 
matplotlib package |127] . 

Appendix A. Notation in HARM and Valencia formulations 

The conservative to primitive inversion routine found in the Einstein Toolkit was 
originally developed as part of the HARM GRMHD code [63j, and reflects its notation 
internally. For historical reasons the development of GRMHD in the Einstein Toolkit 
follows the Valencia formulation of equations P [TTHSUI 175] , and this notation was used 
throughout the first paper summarising the Einstein Toolkit [22] • To convert from one 
set of notation to the other, we note that the equivalence of various physical quantities 
in Table I 



Appendix B. Comparison of reconstruction methods 



We briefly compare the new reconstruction methods (section 4.2) that have been 
implemented in GRHydro. In [3B] a detailed comparison between ePPM and oPPM 
for cell-centred and vertex-centred AMR was performed. These authors found 
that ePPM is superior to oPPM in all cases they considered, and especially in 
those using cell-centred AMR. Here, we compare oPPM, ePPM, WEN05, and MPS 
reconstruction using the example of an magnetised dynamically evolved TOV star 
using vertex-centred AMR. We use 7 levels of mesh refinement with cubical extent 
[4OOM0,2OOM0,5OM0,3OM0,19M0,12M0] and employ octant symmetry. The finest 
level has a resolution of Ax = O.2M0 and contains the entire star. We use the HLLE 
Riemann solver. We evolve spacetime with a F-driver gauge parameter of r; = 0.5 
and 5th-order Kreiss-Oliger dissipation with e = 0.1 is added to the right-hand side of 
spacetime variable. We integrate in time with a fourth-order Runge-Kutta integrator 
and a Courant factor of 0.4. We construct the initial initial data by solving the TOV 
equations using a polytropic EOS with F = 2 and scale parameter K = 100 M0. The 
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Figure Bl. Comparison of oPPM, ePPM, WEN05, and MPS reconstruction based 
on the example of a dynamically evolved TOV star. We show the L2-iiorm of the 
Hamiltonian constraint ||-ff||2- MPS performs best, while oPPM performs worst in 
this simple test. 



central density is set to pc 
EOS witli r = 2. 



1.28 X 10" 



■^o'- 



During tlie evolution, we use a F-law 



In figure Bl , we show the L2-norni of the Hamiltonian constraint | |-ff 1 12 as a function 



of time. All simulation settings are identical, except for the reconstruction method. We 
find that MPS reconstruction performs better than ePPM or WEN05 in this simple 
test. WEN05 performs worse than ePPM and MPS. The worst results are obtained 
with oPPM which reduces to first order at the smooth density maximum at the center 
of the star. In the vertex-centred case, however, the error due to oPPM reconstruction 
is not as large as in the cell-centred case [BB], since in the former, the central density 
maximum is approximately located on a grid point. With oPPM, the dominant source 
of error is thus generated at the center of the star. The other reconstruction methods 
maintain high order at the smooth density maximum. With them, the dominant source 
of error is generated at the steep density drop at the surface of the star. We emphasise, 
however, that this test is only one of many possible tests. A more detailed comparison 
would be necessary for more general conclusions. 



Appendix C. Full set of planar MHD shocktubes 



The five test cases collected in [5B] are widely used to test new MHD codes. We verified 
that our code passes all tests and compare the results to the exact solution of |13] • The 



main text |5.2| contains detailed description of our test setup. 

In figures [UT] and [^2] we present results for the moderate and strong blast wave tests 
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Figure CI. Evolution of the Balsara2 shocktube case, performed with divergence 
cleaning, with all conventions as in figure [3J The numerical code reproduces the 
analytical result very well, with only some oscillations visible in the density p near 
the shock front, capturing the sharp shock features at least as well as |15]. A more 
detailed description of the test setup and parameters can be found in the main text in 
section [ 



of [nS] which prove much more challenging to the code due to the Lorentz contraction 



of the fast moving shock. Finally in figure [C3] we present the non-coplanar test problem 
of 
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Figure C2. Evolution of the BalsaraS shocktube case, performed with divergence 
cleaning, with all conventions as in figure [3] This test proved more demanding on 
our code, in particular the narrow right moving fast shock is not resolved well using 
1600 grid points, resulting in a sizeable overshoot near the shock. This feature is 
common to other codes [SS] independent of the specific formulation used. A more 
detailed description of the test setup and parameters can be found in the main text in 
section 15.21 
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Figure C3. Evolution of the BalsaraS siiocktube case, performed with divergence 
cleaning, with all conventions as in figure |3] A more detailed description of the test 
setup and parameters can be found in the main text in section |5.2[ 
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Appendix D. Monopole tests modelling high frequency noise 

We also include initial data in our tests that models high-frequency numerical noise in 



the magnetic field. Here, the spatial dependence of B^ is set as in section |5.1| and the 
other components of the i?-field are set to the same amplitude or to zero, with non- 
zero components multiplied by (— l)^^*'-?'^)^ where N(i,j, k) is a function of the indices 
i, j, and k of the grid cell in the x, y, and 2;-directions, to produce an alternating 
pattern. For a one- dimensional alternating pattern we set N{i,j,k) = i + j + k, and 
only set B^ to be non-zero and alternating. For a two-dimensional test, we use the same 
function for N{i,j, k) and apply it to B^ to make it alternate, but set B^ to the original 
Gaussian state and B^ = 0. For a three-dimensional test, we set N{i,j,k) = i + j, so 
the alternation is two-dimensional in character, but set B^ = B^ = B^ and allow all 
three to alternate. 



Figure |D1| shows the performance of the divergence cleaning algorithm in removing 
high frequency constraint violating noise such as noise generated during the numerical 
evolution. Clearly the algorithm is very efficient in damping out the high frequency 
constraint violations, reducing them much more quickly than the simple Gaussian pulse 
present in figure [2j 
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Figure Dl. Behaviour of our divergence cleaning scheme, demonstrated by monopole 
damping and advection for a three-dimensional alternating Gaussian configuration. 
Conventions are the same as in figure [2] We still see more significant divergences 
travelling further from their initial location for k = 1, but damping is markedly faster 
for the case k = 100 when dealing with high-frequency noise as compared to the 
lower- frequency modes characterising the Gaussian monopole, seemingly because of 
cancellation effects from neighbouring points. The choice k = 10 still yields the best 
results for numerically eliminating spurious divergence of the magnetic field. 
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